home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / xdiv.cc < prev    next >
C/C++ Source or Header  |  1996-10-11  |  8KB  |  374 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #ifdef HAVE_CONFIG_H
  24. #include <config.h>
  25. #endif
  26.  
  27. #include <cassert>
  28.  
  29. #include "CMatrix.h"
  30. #include "dMatrix.h"
  31. #include "oct-cmplx.h"
  32.  
  33. #include "error.h"
  34. #include "xdiv.h"
  35.  
  36. static inline int
  37. result_ok (int info, double rcond, int warn = 1)
  38. {
  39.   assert (info != -1);
  40.  
  41.   if (info == -2)
  42.     {
  43.       if (warn)
  44.     warning ("matrix singular to machine precision, rcond = %g", rcond);
  45.       else
  46.     error ("matrix singular to machine precision, rcond = %g", rcond);
  47.  
  48.       return 0;
  49.     }
  50.   else
  51.     return 1;
  52. }
  53.  
  54. template <class T1, class T2>
  55. bool
  56. mx_leftdiv_conform (T1 a, T2 b)
  57. {
  58.   int a_nr = a.rows ();
  59.   int b_nr = b.rows ();
  60.  
  61.   if (a_nr != b_nr)
  62.     {
  63.       int a_nc = a.cols ();
  64.       int b_nc = b.cols ();
  65.  
  66.       gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc);
  67.       return false;
  68.     }
  69.  
  70.   return true;
  71. }
  72.  
  73. template bool mx_leftdiv_conform (const Matrix&, const Matrix&);
  74. template bool mx_leftdiv_conform (const Matrix&, const ComplexMatrix&);
  75. template bool mx_leftdiv_conform (const ComplexMatrix&, const ComplexMatrix&);
  76. template bool mx_leftdiv_conform (const ComplexMatrix&, const Matrix&);
  77.  
  78. template <class T1, class T2>
  79. bool
  80. mx_div_conform (T1 a, T2 b)
  81. {
  82.   int a_nc = a.cols ();
  83.   int b_nc = b.cols ();
  84.  
  85.   if (a_nc != b_nc)
  86.     {
  87.       int a_nr = a.rows ();
  88.       int b_nr = b.rows ();
  89.  
  90.       gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
  91.       return false;
  92.     }
  93.  
  94.   return true;
  95. }
  96.  
  97. template bool mx_div_conform (const Matrix&, const Matrix&);
  98. template bool mx_div_conform (const Matrix&, const ComplexMatrix&);
  99. template bool mx_div_conform (const ComplexMatrix&, const ComplexMatrix&);
  100. template bool mx_div_conform (const ComplexMatrix&, const Matrix&);
  101.  
  102. // Right division functions.
  103. //
  104. //       op2 / op1:   m   cm
  105. //            +--   +---+----+
  106. //   matrix         | 1 |  3 |
  107. //                  +---+----+
  108. //   complex_matrix | 2 |  4 |
  109. //                  +---+----+
  110.  
  111. // -*- 1 -*-
  112. Matrix
  113. xdiv (const Matrix& a, const Matrix& b)
  114. {
  115.   if (! mx_div_conform (a, b))
  116.     return Matrix ();
  117.  
  118.   Matrix atmp = a.transpose ();
  119.   Matrix btmp = b.transpose ();
  120.  
  121.   int info;
  122.   if (btmp.rows () == btmp.columns ())
  123.     {
  124.       double rcond = 0.0;
  125.       Matrix result = btmp.solve (atmp, info, rcond);
  126.       if (result_ok (info, rcond))
  127.     return Matrix (result.transpose ());
  128.     }
  129.  
  130.   int rank;
  131.   Matrix result = btmp.lssolve (atmp, info, rank);
  132.  
  133.   return result.transpose ();
  134. }
  135.  
  136. // -*- 2 -*-
  137. ComplexMatrix
  138. xdiv (const Matrix& a, const ComplexMatrix& b)
  139. {
  140.   if (! mx_div_conform (a, b))
  141.     return ComplexMatrix ();
  142.  
  143.   Matrix atmp = a.transpose ();
  144.   ComplexMatrix btmp = b.hermitian ();
  145.  
  146.   int info;
  147.   if (btmp.rows () == btmp.columns ())
  148.     {
  149.       double rcond = 0.0;
  150.       ComplexMatrix result = btmp.solve (atmp, info, rcond);
  151.       if (result_ok (info, rcond))
  152.     return result.hermitian ();
  153.     }
  154.  
  155.   int rank;
  156.   ComplexMatrix result = btmp.lssolve (atmp, info, rank);
  157.  
  158.   return result.hermitian ();
  159. }
  160.  
  161. // -*- 3 -*-
  162. ComplexMatrix
  163. xdiv (const ComplexMatrix& a, const Matrix& b)
  164. {
  165.   if (! mx_div_conform (a, b))
  166.     return ComplexMatrix ();
  167.  
  168.   ComplexMatrix atmp = a.hermitian ();
  169.   Matrix btmp = b.transpose ();
  170.  
  171.   int info;
  172.   if (btmp.rows () == btmp.columns ())
  173.     {
  174.       double rcond = 0.0;
  175.       ComplexMatrix result = btmp.solve (atmp, info, rcond);
  176.       if (result_ok (info, rcond))
  177.     return result.hermitian ();
  178.     }
  179.  
  180.   int rank;
  181.   ComplexMatrix result = btmp.lssolve (atmp, info, rank);
  182.  
  183.   return result.hermitian ();
  184. }
  185.  
  186. // -*- 4 -*-
  187. ComplexMatrix
  188. xdiv (const ComplexMatrix& a, const ComplexMatrix& b)
  189. {
  190.   if (! mx_div_conform (a, b))
  191.     return ComplexMatrix ();
  192.  
  193.   ComplexMatrix atmp = a.hermitian ();
  194.   ComplexMatrix btmp = b.hermitian ();
  195.  
  196.   int info;
  197.   if (btmp.rows () == btmp.columns ())
  198.     {
  199.       double rcond = 0.0;
  200.       ComplexMatrix result = btmp.solve (atmp, info, rcond);
  201.       if (result_ok (info, rcond))
  202.     return result.hermitian ();
  203.     }
  204.  
  205.   int rank;
  206.   ComplexMatrix result = btmp.lssolve (atmp, info, rank);
  207.  
  208.   return result.hermitian ();
  209. }
  210.  
  211. // Funny element by element division operations.
  212. //
  213. //       op2 \ op1:   s   cs
  214. //            +--   +---+----+
  215. //   matrix         | 1 |  3 |
  216. //                  +---+----+
  217. //   complex_matrix | 2 |  4 |
  218. //                  +---+----+
  219.  
  220. Matrix
  221. x_el_div (double a, const Matrix& b)
  222. {
  223.   int nr = b.rows ();
  224.   int nc = b.columns ();
  225.  
  226.   Matrix result (nr, nc);
  227.  
  228.   for (int j = 0; j < nc; j++)
  229.     for (int i = 0; i < nr; i++)
  230.       result (i, j) = a / b (i, j);
  231.  
  232.   return result;
  233. }
  234.  
  235. ComplexMatrix
  236. x_el_div (double a, const ComplexMatrix& b)
  237. {
  238.   int nr = b.rows ();
  239.   int nc = b.columns ();
  240.  
  241.   ComplexMatrix result (nr, nc);
  242.  
  243.   for (int j = 0; j < nc; j++)
  244.     for (int i = 0; i < nr; i++)
  245.       result (i, j) = a / b (i, j);
  246.  
  247.   return result;
  248. }
  249.  
  250. ComplexMatrix
  251. x_el_div (const Complex a, const Matrix& b)
  252. {
  253.   int nr = b.rows ();
  254.   int nc = b.columns ();
  255.  
  256.   ComplexMatrix result (nr, nc);
  257.  
  258.   for (int j = 0; j < nc; j++)
  259.     for (int i = 0; i < nr; i++)
  260.       result (i, j) = a / b (i, j);
  261.  
  262.   return result;
  263. }
  264.  
  265. ComplexMatrix
  266. x_el_div (const Complex a, const ComplexMatrix& b)
  267. {
  268.   int nr = b.rows ();
  269.   int nc = b.columns ();
  270.  
  271.   ComplexMatrix result (nr, nc);
  272.  
  273.   for (int j = 0; j < nc; j++)
  274.     for (int i = 0; i < nr; i++)
  275.       result (i, j) = a / b (i, j);
  276.  
  277.   return result;
  278. }
  279.  
  280. // Left division functions.
  281. //
  282. //       op2 \ op1:   m   cm
  283. //            +--   +---+----+
  284. //   matrix         | 1 |  3 |
  285. //                  +---+----+
  286. //   complex_matrix | 2 |  4 |
  287. //                  +---+----+
  288.  
  289. // -*- 1 -*-
  290. Matrix
  291. xleftdiv (const Matrix& a, const Matrix& b)
  292. {
  293.   if (! mx_leftdiv_conform (a, b))
  294.     return Matrix ();
  295.  
  296.   int info;
  297.   if (a.rows () == a.columns ())
  298.     {
  299.       double rcond = 0.0;
  300.       Matrix result = a.solve (b, info, rcond);
  301.       if (result_ok (info, rcond))
  302.     return result;
  303.     }
  304.  
  305.   int rank;
  306.   return a.lssolve (b, info, rank);
  307. }
  308.  
  309. // -*- 2 -*-
  310. ComplexMatrix
  311. xleftdiv (const Matrix& a, const ComplexMatrix& b)
  312. {
  313.   if (! mx_leftdiv_conform (a, b))
  314.     return ComplexMatrix ();
  315.  
  316.   int info;
  317.   if (a.rows () == a.columns ())
  318.     {
  319.       double rcond = 0.0;
  320.       ComplexMatrix result = a.solve (b, info, rcond);
  321.       if (result_ok (info, rcond))
  322.     return result;
  323.     }
  324.  
  325.   int rank;
  326.   return a.lssolve (b, info, rank);
  327. }
  328.  
  329. // -*- 3 -*-
  330. ComplexMatrix
  331. xleftdiv (const ComplexMatrix& a, const Matrix& b)
  332. {
  333.   if (! mx_leftdiv_conform (a, b))
  334.     return ComplexMatrix ();
  335.  
  336.   int info;
  337.   if (a.rows () == a.columns ())
  338.     {
  339.       double rcond = 0.0;
  340.       ComplexMatrix result = a.solve (b, info, rcond);
  341.       if (result_ok (info, rcond))
  342.     return result;
  343.     }
  344.  
  345.   int rank;
  346.   return a.lssolve (b, info, rank);
  347. }
  348.  
  349. // -*- 4 -*-
  350. ComplexMatrix
  351. xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b)
  352. {
  353.   if (! mx_leftdiv_conform (a, b))
  354.     return ComplexMatrix ();
  355.  
  356.   int info;
  357.   if (a.rows () == a.columns ())
  358.     {
  359.       double rcond = 0.0;
  360.       ComplexMatrix result = a.solve (b, info, rcond);
  361.       if (result_ok (info, rcond))
  362.     return result;
  363.     }
  364.  
  365.   int rank;
  366.   return a.lssolve (b, info, rank);
  367. }
  368.  
  369. /*
  370. ;;; Local Variables: ***
  371. ;;; mode: C++ ***
  372. ;;; End: ***
  373. */
  374.